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ABSTRACT 

We describe an algorithm to generate temperature and polarization maps 
of the cosmic microwave background radiation containing non-Gaussianity of 
arbitrary local type. We apply an optimized quadrature scheme that allows 
us to predict and control integration accuracy, speed up the calculations, and 
reduce memory consumption by an order of magnitude. We generate 1000 non- 
Gaussian CMB temperature and polarization maps up to a multipole moment of 
^max = 1024. We validate the method and code using the power spectrum and 
the fast cubic (bispectrum) estimator and find consistent results. The simulations 
are provided to the communitjj^j 
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1. Introduction 

The simplest models of inflation predict almost perfectly Gaussian primordial fluctu- 



ations, generated by a single scalar quantum field in ground state (Guth 1981 Bardeen 



et al. 1983 Mukhanov et al. 1992), but a large number of alternative scenarios can easily be 



constructed. To test competing inflationary models, measurements of statistical properties 
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of the cosmic microwave background (CMB) radiation have turned out to be of particular 
importance. Combined with constraints on the scalar spectral index n s and the search for 
gravitational waves imprinted on the polarization signature, the test for non-Gaussianity is 
a fundamental means to probe the physical processes of inflation. 

Among all inflationary models predicting significant levels of non-Gaussianity, two broad 
classes can be distinguished. Non-Gaussianity of equilateral type is realized primarily in 



models with non-minimal Lagrangian including higher order derivatives (Alishahiha et al. 



2004 Senatore 2005; Chen 2005 Langlois et al. 2008). Non-Gaussianity of local type is 



achieved to very good approximation in multi-field inflation (Moroi & Takahashi 2001 En- 



qvist & Sloth 


2002 


Lyth et al. 


2003 


), or in cyclic/ekpyrotic 


2001 


Steinhardt h Turok 


2002 


Lehners & Steinhardt 2008) 



Concentrating on local non-Gaussianity, we parameterize the primordial curvature per- 
turbations, $, by introducing an additional quadratic dependence on a purely Gaussian 



auxiliary field $l, that is local in real space, of the form (Verde et al. 2000 Komatsu & 

$(r) = $ L (r) + /nl $nl(>) , 



Spergel 2001) 



where $nl(j) is defined as 



$ NL (r) = $£(r) 



1) 



(2) 



and /nl is the dimensionless measure of the amplitude of non-Gaussianity. Primordial non- 
Gaussianity in the curvature perturbations $ will be encoded in the CMB signal. 

Simulations of maps containing non-Gaussianity of local type have been extensively used 
in the context of WMAP data analysis. They play a crucial role in probing the sensitivity 
of analysis tools and provide the opportunity to quantitatively estimate the contribution 
of secondary anisotropies or instrumental effects to the measured level of non-Gaussianity 
in experimental data. With the advent of Planck, probing the non-Gaussian contribution 
within the CMB radiation even more accurately, the requirements for high resolution, high 
accuracy simulations of non-Gaussian CMB temperature and polarization maps will further 
increase. 

To meet the demand for simulated non-Gaussian maps, several different approaches 
have been taken. First simulations of temperature maps with primordial non-Gaussianity of 



local type have generated the underlying primordial perturbation in Fourier space (Komatsu 



et al.|200~3 ). This approach is computationally very demanding while it is difficult to preserve 



numerical accuracy. A different method has been proposed in Liguori et al. (2003, 2007), 
where the authors work with 'filter' functions to introduce the proper spatial correlations 
of the primordial potential. Recently, a fast, specifically tailored algorithm for the weakly 
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non-Gaussian regime has been introduced by Smith & Zaldarriaga (2006), that focuses 
on simulating maps with a given three-point function. While it is not restricted to non- 
Gaussianity of local type, higher order correlations are not guaranteed to match the model. 



The algorithm presented here was closely inspired by the work of Liguori et al. (2003). 
We focus on an enhancement of their algorithm in view of its numerical efficiency Our idea is 
to precompute quadrature nodes and weights; this is similar in spirit to|Smith Zaldarriaga 



(2006), but aims at assuring accurately simulated maps to all correlation orders, rather than 



focusing exclusively on the three-point function. 

This paper is organized as follows. In Sect. [2j we present a new approach to simulate 
non-Gaussian temperature and polarization maps. An optimization scheme is provided in 
Sect. [3] that allows for an increase in computational efficiency We then apply the fast 
estimator to simulated CMB maps to check our results for consistency (Sect. [1]). Finally, we 
summarize our findings in Sect. |5j 

Throughout the paper we assume the following M^MAP5+BAO+SN cosmological pa- 
rameter^ = 0.721, Q c h 2 = 0.1143, Q b h 2 = 0.02256, A^(0.002 Mpc' 1 ) = 2.457- 10~ 9 , 
h = 0.701, n s = 0.96, and r = 0.084. 



2. Simulation of non— Gaussian CMB maps 

We describe a new, direct method to simulate non-Gaussian CMB temperature and 
polarization maps below. Our objective is to generate a set of linear and non-linear spherical 
harmonic coefficients that are valid realizations of temperature and polarization fluctuations, 
{ol im, ^nl em}, for a given cosmological model. A map with any desired level of non- 
Gaussianity, /nl, can then be realized by linear combination, 

&lm = «L £m + /nl ' a NL im ■ (3) 

The expansion coefficients ae m of the CMB temperature and polarization anisotropies in 



harmonic space are related to the primordial fluctuations &£ m (k) via the equation (Komatsu 



et al. 2003) 



«e m ~ -i^r I dk k 2 $im{k) g}(k) . (4) 

Here, g\{k) is the transfer function of temperature (i = T) or polarization (i = E) in mo- 
mentum space. Analogously, we can define an equivalent equation as a function of comoving 
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distance, 

drr 2 $tm(r) a\(r) , (5) 



where we have used the real space transfer function according to 

4(r) = l J dkk 2 gi(k)j e (kr), (6) 

where je{kr) denotes the spherical Bessel function of order i. 

We can now outline our recipe for simulating non-Gaussian CMB maps as following: 
(i) Generate the multipole moments of a purely Gaussian gravitational potential $l tmif) as 
a function of conformal distance, (ii) Compute the spherical harmonic transform to derive 
the corresponding expression in pixel space, ^\jr). (iii) Square it and subtract the variance 
according to Eq. [2] to get the non-Gaussian potential $ NL (r). (iv) Inverse transform to 
spherical harmonic space to obtain $ NL e m ( r )- ( v ) Solve the radial integral Eq. [5]for $ L tm{ r ) 
and $nl tm(r) separately to compute {a£ im , af im , a£ L im , af L im }. 

One difficulty in this approach is that we have to take into account the radial corre- 
lation of the gravitational potential in step (i). Its covariance matrix is determined by the 
primordial power spectrum predicted by inflation, V(k), and is given by ( |Liguori et aLp003] ) 



f A 2 (k) 

<$l < imi (n) $l i 2 m 2 (r2)) = 4vr 5% S% J dk -f^- j^kv,) j h (kr 2 ) , (7) 
where we have replaced V{k) by 

^Uk) = ^-V(k), (8) 



2tt 2 

that is constant for vanishing spectral tilt (n s = 1). The covariance matrix will be denoted 
by -P$ e{ r i, r 2) in what follows. To draw a random realization of the linear gravitational 
potential at distances r = (r 1; r 2 , . . . , r n ), we calculate 

$Le m (r) = Pl /2 c9, (9) 

where g is a vector of independent complex Gaussian random variables with zero mean and 
unit variance. 

For this algorithm to run efficiently, we have to reduce the number of quadrature points 
in the numerical evaluation of the radial integral (Eq. |5|, to keep the number of computation- 
ally expensive spherical harmonic transformations necessary to generate the non-Gaussian 
gravitational potential as low as possible. Details of the implementation together with an 
optimization scheme will be described in the next section. 
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3. Implementation and Optimization 

To be able to perform the steps outlined in the last section, we first have to precompute 
the necessary auxiliary data. This needs to be done only once for a given set of cosmological 
parameters. First, we obtained the transfer functions in momentum space from a modified 
version of the latest CAMB software packagd^ (Lewis et al. 2000). We then derived their 



equivalent expressions in real space using Eq. |6j Examples of temperature and polarization 
transfer functions as a function of conformal distance for several multipole moments are 
shown in Fig. [3] 

As a next step, we calculate the covariance matrix of the gravitational potential on a 
fine grid with N to t = 400 shells from the origin to the present time cosmic horizon (Eq. [7]). 
As a start, we resolve the last scattering surface with uniform spacing using an increment 
of Ar « 3.5 Mpc and chose a larger interval elsewhere (Ar w 100 Mpc). This simple 
approach will be refined later. Using the derived quantities, it is now possible to generate 
$l tm{ r i)i ^nl im{ r i) an d numerically solve the radial integral Eq. [5] to obtain simulated 
non-Gaussian CMB maps. However, significant improvement in the numerical evaluation 
of the integral is achievable by choosing both weights and quadrature points in an optimal 
way, as we will show in the following. 

Keeping the multipole moment (£, m) fixed for simplicity, we want to accurately compute 
the integral 

1= drr 2 a(r) $(r) . (10) 



This is done in a two step process: We first approximate Eq. 10 with a discrete sum over N tot 
shells. Then, we try to obtain comparable accuracy with fewer shells N <C N tot introducing 
weights, 

N 

/ = $(r 4 ) , (11) 

i=l 

where the gravitational potential is evaluated at the nodes r^ and weighted by the factors 
W{. Now, we can derive the expectation value of the quadratic error 

((/-/) 2 > = X; A, ^Jdrr 2 a(r)Mr)J -w T P*w, (12) 

where we have introduced the eigenvalues and eigenvectors <pk of the covariance matrix 
of the potential on the fine grid with Ntot — 400 elements. We show for several multipole 



3 Obtaincd from http://camb.info 
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mo merits in the left panel of Fig. [3] If the eigenvalues decrease sufficiently fast, the error is 
expected to be low already for a small number of quadrature points N. This seems especially 
to be true on large angular scales. However, this finding is partially counterbalanced by 
the fact that the transfer functions are significantly different from zero at small radii for 
low multipole moments (late ISW effect, reionization) , enforcing the inclusion of additional 
nodes. 



Based on the expression for the expected quadratic error, it is straightforward to calcu- 

a 

t)ir, 



late optimal weights by satisfying the condition jf-((I — I) 2 ) = 0, which leads to a system 



of iV linear equations, 

N N t . 



1 * J *tOt /» 

Y] p *jj(n,rj) Wj = ^2 h 4>k{ri) / drr 2 a(r) <p k {r) 
j=l k=i -> 



(13) 



Even more important, Eq. 12 allows us to formulate a greedy algorithm to compute 
optimal quadrature points. We select a subset of nodes out of the fine radial grid with 
400 elements iteratively, in each step including the point that most efficiently reduces the 
remaining error. To simultaneously optimize for temperature and polarization, we add the 
expectation values of the two errors with equal weights. We use the outcome of the procedure 
to tune the radii of the input grid with 400 elements. We choose a smaller spacing down 
to Ar = 1.2 Mpc at the last scattering surface, where nodes were selected with the highest 
priority, and a larger step size up to Ar = 140 Mpc at distances, where the quadrature points 
were classified as less important. Then, we repeated the optimization process a second time. 
In the right panel of Fig. [3j we visualize the first 100 iterations of the optimization scheme. 
We display the expectation value of the relative quadratic error for N = 30, 50 and N = 70 
quadrature points in Fig. [3j The raise in error towards the largest angular scales is caused 
by the increasing contribution from late ISW effect and reionization. 

As a last step, we show how to reduce the memory consumption of the code. Whereas the 
potentials $l e m ( r ) and $nl £m( r ) can be co-added to the spherical harmonic coefficients 
of the CMB map a\\ lm and a l NL tm (Eq. [5]) one shell after another, the radial correlation 
of < &l(^) forces one to generate it at all N shells simultaneously, potentially requiring large 
amounts of memory. To circumvent this problem, we keep the random seeds that were used to 
draw the potential. By means of the seeds, we are able to easily regenerate the gravitational 
potential at any radii (fx, r 2 , . . . , rjy). Thus, we only store its real space representation at the 
radius that is currently added to the CMB map, substantially reducing the overall memory 
consumption of the algorithm. 

Having optimized the simulation algorithm in this way, we generated N sim = 1000 
realizations of temperature and polarization CMB maps. We chose a HEALPix resolution 
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Fig. 1. — Real space transfer functions. We show examples of the real space transfer func- 
tions of temperature (left panel) and polarization (right panel) for three different multipole 
moments £ = 5, 20, and 100. At low £, the effect of the late time ISW effect is clearly visible. 
Reionization occurred at about r = 10 Gpc. 




Fig. 2. — Optimization scheme. Left panel: We display the largest 100 eigenvalues of the 
covariance matrix of $l e(f), normalized and in descending order for i = 5, 20, and 100. For 
low multipole moments, the number of quadrature points can be reduced most efficiently 
Right panel: The radial positions of the shells included in the first 100 iteration steps. For 
illustrative purposes, we interchanged open and filled symbols every 10 iterations. The most 
important nodes are included first. 
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parameter of n si( i e = 512 and a maximum multipole moment of £ ma x — 1024. We used 
N = 70 quadrature points for evaluation of Eq. [5j although we stress that this choice is 
conservative and it is possible to derive reasonable results with smaller values of N. 

With these input parameters, we aim for sub-percentage accuracy of the final map over 
the entire range of multipole moments, guided by the intrinsic precision of the underlying 
transfer functions, running CAMB with RECFAST at standard accuracy An example is 
displayed in Fig. |3j where we illustrate a realization of temperature and polarization maps 
of the linear and non-linear part of the CMB. We show the averaged power spectra of all 
simulations along with a comparison to the theoretical values in Fig. [3j A detailed compar- 
ison to the expected statistical fluctuations (oc -/V sim ' ) reveals remaining slight systematic 
deviations for the TT and EE spectra at high i at the level of less than 1% of the input 
power spectrum. If required this error could be further reduced by adding integration nodes. 
It takes about 20 minutes to generate a single map with the given resolution on a single Intel 
Xenon processor with a clock rate of 2.33 GHz, requiring only a modest amount of memory 
(« 400 MB). The most time consuming part is the evaluation of the spherical harmonic 
transforms necessary to compute the non-Gaussian potential. 

The algorithm described here generates valid realizations of primordial curvature per- 
turbations in real space. This itself is an interesting quantity and can be used to e.g. test the 
performance of reconstruction techniques as we will show in the next section. We visualize 
the 3D gravitational potential $L( r ) and $NL( r ) in Fig. [3] long-distance correlations on 
large scales are in evidence. 



Our simulation algorithm is conceptually very similar to the method proposed in Liguori 



et al. (2007), where the authors generate the gravitational potential on 400 shells, requiring 
800 spherical harmonic transforms to calculate a single non-Gaussian CMB map, and report 
a runtime of 3 hours for £ max = 500. By applying our optimized quadrature scheme, we have 
demonstrated that it is possible to reduce the number of transforms considerably, resulting in 
an increase of computational efficiency Another, albeit more formal difference is the way the 
gravitational potential is generated. We use the real space covariance matrix to draw 



directly, whereas the authors of Liguori et al. (2007) compute the gravitational potential by 



performing an integral over uncorrelated random numbers weighted by 'filter' functions. 



In Smith & Zaldarriaga (2006), where the authors focused on a perturbative reproduc- 
tion of the correct bispectrum in the regime of weak non-Gaussianity, a runtime of about 
three minutes is reported to simulate one non-Gaussian CMB temperature map at an an- 
gular resolution of £ ma x = 1000. Although slower by an order of magnitude, and tuned for 
local non-Gaussianity, the algorithm presented here is capable of simulating both temper- 
ature and polarization maps (i.e. three maps for the stokes parameter I, Q, and U) within 
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the same framework and with nearly the same computational cost compared to temperature 
alone. Furthermore, as recently pointed out by Hanson et al. (2009), in the case of local 



non-Gaussianity an additional modification of the algorithm of Smith & Zaldarriaga (2006) 



is necessary to suppress the power spectrum of the non-Gaussian part of a simulated map, 
found to be artificially enhanced by several orders of magnitudes on large angular scales. 

Notwithstanding the aforementioned higher computational costs, we regard our method 
as useful for the study of local non-Gaussianity, because the simulated maps are well suited 
to test any kind of estimator, e.g. based on Minkowsky functionals (Spergel et al. 2007 



Hikage et al. 2008), or a wavelet analysis (Martinez-Gonzalez et al. 2002 Mukherjee & 



Wang 2004). If a detection of nonzero /nl is reported, it will be important to confirm the 



result with alternative statistical tools, as they are sensitive to different systematic effects. 



In the following section, we apply the KSW estimator (Komatsu et al. 2005) to our set 



of simulated maps with known non-Gaussian contribution to test whether the input values 
for / NL can be recovered. 



4. Bispectrum Analysis 



As we do not aim to describe the fast estimator in detail, we include a brief summary 
here and refer the reader to the extensive literature for further details and a comprehensive 



Komatsu et al. 


2005 


Smith & Zaldarriaga 


2006 


Creminelli et al. 


2007 



Yadav et al.||2007| ). 



To estimate the non-Gaussianity of a CMB map, one constructs the statistic S pr i m out 
of a cubic combination of the data, 



S 



prim 



drr 2 / d 2 h A(r,n) B 2 (r,n) 



The radial integral runs over two filtered maps, 

A(r,n)= £ £(0?«S(r) 4m YtUn) 

i,j=T,E e,m 



i,j=T,E l,m 

that are constructed using the auxiliary functions 

2 

7T 



a at] 



dk k g\(k) je(kr) 



/3j( r ) = 4tt / dk 



k 



g\{k) je(kr) 



(14) 

(15) 
(16) 

(17) 
(18) 
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and the inverse of the matrix containing the CMB power spectrum elements, 

-l 

nTE nEE 



C7 1 



(19) 



One of these maps, B(r,n), is exactly the Wiener filter reconstruction of the underlying 
gravitational potential $(r). With the simulation algorithm presented in Sect. [2j it is pos- 
sible to compare the potential used to synthesize the map with its reconstruction directly. 
An example is shown in Fig. [4], where we depict the reconstruction of the potential around 
last scattering using only temperature, and using both, temperature and polarization infor- 
mation. 

As the estimator S pr i m is proportional to the non-Gaussianity parameter we can 
calculate its expectation value by applying a suitable normalization, 



NL 



i,j,k,o,p,q h<(-2<h '1*2*3 
=T,E 



(C-% {C-'t (C" J )J B, 



prim ; 



where = 6, when l\ — £ 2 = £3, 2, when £\= £2^ £3 or £\ 7^ £ 2 

We further introduced the theoretical bispectrum for / NL 



(20) 

?3, and 1 otherwise. 



1, B^'J* 1 " 1 , which is defined as 



Kiir = 2 4« a J drr' [^(r)^ 2 (r)a|(r) + ^ 3 (r)^(r)^ 2 (r) + /3i 2 (r)^ 3 (r)^(r)] , 

(21) 

where the prefactor is given by 



(2£ 1 + 1)(2£ 2 + 1)(2£ 3 + 1) / £1 £2 h 
4tt V 



(22) 



We used the equations above to implement the fast estimator for temperature and 
polarization. As our primary goal is to validate our simulation algorithm, we do not take 
into account possible instrumental effects, sky cut, or noise. To test our simulations, we 
generate two sets of 1000 CMB temperature and polarization maps with resolution parameter 
n S ide = 512, and £ max = 1024. We consider one sample of purely Gaussian realizations of 
the CMB sky (/ NL = 0), and one non-Gaussian sample with a fiducial value of /nl = 100. 
We then run the fast estimator on the maps to compute an estimate of /nl- We show the 
distribution of the derived values in Fig. |4j We find the input parameters to be recovered 
well, the means of the distributions are (f§ L ) = —0.1 and (/nlT) = 98.4 for the Gaussian and 
non-Gaussian simulations, respectively. The estimated standard deviations are <Tjg^ =2.4 
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and a 



■i N 



8.4, compared to the expected error predicted from a Fisher information matrix 



analysis of <3 F j^ eT = 2.4. We conclude that the algorithm outlined in Sect. 2]and implemented 
as described in Sect. [3] produces valid realizations of non-Gaussian CMB temperature and 
polarization maps. 



5. Summary 

In this paper, we introduced a new algorithm to simulate temperature and polarization 
CMB maps containing non-Gaussianity of arbitrary local type. In the proposed scheme, 
we generate spherical harmonic coefficients of the Gaussian potential as a function of con- 
formal distance, taking into account the proper radial correlations. Then, the potential is 
transformed to pixel space to compute the associated non-Gaussian contribution. Finally, 
we make use of the real space representation of the transfer functions to perform the line 
of sight integral in order to calculate Gaussian and non-Gaussian contribution to the CMB 
maps. 

We developed and applied a quadrature scheme that allows us to increase the numerical 
efficiency of the code. As a starting point, we derived an expression to quantitatively calculate 
the mean error introduced by replacing the radial integral by a finite sum. On that basis, 
we were able to choose both, nodes and weights for numerical quadrature in an optimal way. 
As a last step, we successfully reduced the memory consumption of the algorithm. 

For 1/KMAP5+BAO+SN cosmological parameters, we simulated 1000 realizations of 
non-Gaussian CMB temperature and polarization maps with resolution parameters n S id e = 
512 and £ max = 1024. To validate the algorithm, we applied the well studied and widely 
accepted fast cubic (bispectrum) estimator to the simulations. For both, a set of Gaussian 
and non-Gaussian realizations of CMB sky maps, the input parameters were consistently 



recovered. We make our simulations publicly available at http://planck.mpa-garching. 
mpg . de/ cmb/fnl- simulations 



Some of the results in this paper have been derived using the HEALPix (Gorski et al. 



2005] ) package. BDW is partially supported by NSF grants AST 0507676 and AST 07-08849. 
BDW gratefully acknowledges the Alexander v. Humboldt Foundation's Friedrich Wilhelm 
Bessel Award. BDW thanks the Caltech Astrophysics group for their hospitality while this 
work was being completed. 
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Fig. 3. — Error of integration. We depict the relative mean quadratic error introduced by 
approximating the integral Eq. [5]by a sum over N = 30, 50, and 70 elements for temperature 
(left panel) and polarization (right panel). 
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Fig. 4. — Example of simulated non-Gaussian map. We show the linear (left column) and 
the associated non-linear part (right column) of a realization of temperature (first row) 
and polarization intensity (second row) CMB data. The polarization intensity is defined as 
/ = \lQ 2 + U 2 , where Q and U are the Stokes parameters. 
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i 



Fig. 5. — Averaged power spectra. Left panel: We display the power spectra C TT , C EE , 
and C TE of the linear part of the simulated CMB maps, averaged over 1000 simulations. 
We do not show the input power spectra here, as the lines cannot be discerned in this 
view. Right panels: The ratio of the power spectra divided by their theoretical values for 
temperature (XX = TT, upper sub-panel), polarization (XX = EE, middle sub-panel), 
and cross-power spectrum (XX = TE, lower sub-panel) . Oscillatory features in the latter 
are caused by roots of the denominator. The grayish area indicates the 2-a bounds of an 
ideal simulation code. Sub-percentage, systematic deviations for the TT and EE spectra 
remain but are consistent with the precision goal. 
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Fig. 6. — Examples for simulated curvature perturbations. Left panel: We visualize the 
linear gravitational potentials ^^(r), generated on N = 70 shells from the origin (center) 
to the last scattering surface (outermost shells). Right panel: The associated non-Gaussian 
potential, displayed at nonlinear scale. 




Fig. 7. — Wiener filter reconstruction of the gravitational potential. We illustrate the input 
gravitational potentials at the last scattering surface r = 14.0 Gpc (middle panel) 

used to generate a simulated CMB map and its Wiener filter reconstruction based solely on 
temperature data (left panel), and based on both, temperature and polarization data (right 
panel) of the same map. Each patch is 50° on the side. 
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Fig. 8. — Histogram of the recovered / NL values. We display the distribution of estimated / NL 
values when applying the fast estimator to 1000 realization of temperature and polarization 
CMB maps. The input values used for the simulations were /nl = (left panel), and 
/nl = 100 (right panel). 



